Abstract 



Quantum optimal control experiments and simulations have success- 
fully manipulated the dynamics of systems ranging from atoms to biomolecules. 
Surprisingly, these collective works indicate that the effort (i.e., the num- 
ber of algorithmic iterations) required to find an optimal control field 
appears to be essentially invariant to the complexity of the system. The 
present work explores this matter in a series of systematic optimizations of 
the state-to-state transition probability on model quantum systems with 
the number of states A'^ ranging from 5 through 100. The optimizations 
occur over a landscape defined by the transition probability as a function 
of the control field. Previous theoretical studies on the topology of quan- 
tum control landscapes established that they should be free of sub-optimal 
traps under reasonable physical conditions. The simulations in this work 
include nearly 5000 individual optimization test cases, all of which con- 
firm this prediction by fully achieving optimal population transfer of at 
least 99.9% upon careful attention to numerical procedures to ensure that 
the controls are free of constraints. Collectively, the simulation results 
additionally show invariance of required search effort to system dimen- 
sion A''. This behavior is rationalized in terms of the structural features 
of the underlying control landscape. The very attractive observed scaling 
with system complexity may be understood by considering the distance 
traveled on the control landscape during a search and the magnitude of 
the control landscape slope. Exceptions to this favorable scaling behavior 
can arise when the initial control field fiuence is too large or when the 
target final state recedes from the initial state as A'^ increases. 



1 



Exploring Quantum Control Landscapes: 
Topology, Features, and Optimization Scaling 

Katharine W. Moore and Herschel Rabitz 
Department of Chemistry, Princeton University, Princeton, NJ, 08544 

June 10, 2011 



1 Introduction 

The control of quantum phenomena with external fields using optimal control 
theory (OCT) [l|, 0| and optimal control experiments (OCE) Q is currently an 
active area of research [J,|5|. OCT simulations have successfully controlled a va- 
riety of objectives, including state preparation 0,0,01, molecular isomerization 
[sl-fl^. dissociation and orientation/alignment [l7l - [l9} . OCE using ul- 

trafast tailored la ser p ulses have achieved control over many processes including 
state preparation U 21|, selective molecular dissociation [22-24[, generation of 
high order optical harmonics 25l- 27|. and energy transfer and isomerization in 



large biomolecules |28l - l3Cll |. Simulation models consider from 2 to 10 or more 
states, and the atoms/molecules used in OCE often have much larger numbers 
of accessible states. Remarkably, controlling complex quantum systems appears 
to be no more difficult than controlling simple ones, both in simulations and 
experiments, where the level of difficulty is expressed in terms of the number of 
iterations required to converge on the target objective. 

The success of these and other studies suggests that quantum control is gen- 
erally amenable to "easy" solution by optimal search. Recently, the quantum 
control landscape concept was introduced to help rationalize the observed wide 
success of quantum control studies fsi'], where the landscape is defined as the 
functional relationship between the physical objective (e.g., population transfer 
probability Pi^f) and the external control field e{t). Considering a controllable 
target system under reasonable physical assumptions [si*], the topology of the 
dynamical quantum control landscape can be shown to have no suboptimal lo- 
cal maxima or traps 31, Exceptions to this favorable topology have 
been found under unusual circumstances, e.g., when constant control fields e{t) 



are employed (BQ. An important objectfve is to either afiirm the attractive 
theoretical landscape findings or identify the likelihood of encountering land- 
scape traps in the course of typical optimizations under reasonable physical 
conditions. The extensive prior optimal control literature is supportive of the 
landscape theory with often high reported yields [iHl, 0- 1^, 3§- 53 . Such stud- 
ies, however, cannot rigorously assess the landscape topology due to constraints 
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of various types (e.g., control field fluence) limiting access to the highest yields 
on the landscape. Additionally, great numerical care is needed when testing the 
landscape for traps as significant numerical limitations (e.g., insufficient tem- 
poral discretization of the control field) can introduce artificial traps. Thus, in 
the present work we execute a large number of carefully performed numerical 
simulations to assess the ability to climb the landscape without encountering 
traps. 

This work will consider the control objective of maximizing the probability 
Pi^f of population transfer from an initial pure state \i) to some target pure 
state I/) of a closed quantum system undergoing unitary evolution. Although 
in the laboratory the circumstances will typically include additional factors be- 
yond this idealized situation, the objective of maximizing the population in the 
product state is often the ultimate goal. The control objective is to identify a 
suitable field e{t) that maximizes Pi^f at some target time T, which may be 
finite or asymptotic with T oo. Typically, an optimal field is found using 
a suitable search algorithm (see, for example, [1, [lol) to traverse the relevant 
control landscape, which is specified by Pi^f as a functional of the control field, 
Pi_j>/ = Pi^f[e{t)]. Both the global topology and local structure of the control 
landscape may influence the character and duration of the search trajectory 
from an initial (often random) control field to an optimal solution. 

The search effort required to find an optimal control field is an important 
issue for determining the feasibility of performing both quantum control simu- 
lations and experiments, as computational and experimental resources are in- 
evitably limited. In particular, if the effort rises with system complexity, search- 
ing for an optimal control field may become too expensive for complex quantum 
systems. In this work, the complexity of the system is measured by the Hilbert 
space dimension iV, i.e., the number of accessible energy levels of iJo- A large 
body of results from the OCT literature 0, 0, i, [H, [H, [ll [s^jHll performed on 
systems where N ranges from 2 to ~ 10^ suggest that the search effort required 
for population transfer does not scale strongly with N. Although the required 
effort will depend on the convergence criteria, the number of reported algorith- 
mic iterations to achieve convergence is observed to be typically no more than 
^ 10'^, and often ^ 100 or fewer, regardless of N or the particular search algo- 
rithm employed. The invariance of required search effort with respect to TV has 
been numerically demonstrated for the Pi^f objective using so-called kinematic 
control variables (i.e., the elements of the governing unitary transformation, or 
equivalent variables) [6l[. In the present work, the scaling of the effort with N 
to find a solution is systematically studied using dynamic control variables (i.e., 
the control field e{t)) for simple model quantum systems. Here, effort is defined 
as the number of algorithmic iterations required to reach a particular threshold 
value of Pi^f] we put aside the effort per iteration to solve the Schrodinger 
equation, which is strongly dependent on N. This condition corresponds to 
the laboratory situation, where the effort of performing an experiment is not 
necessarily dependent on the complexity of the target molecule. The dynami- 
cal control findings throughout the paper will be compared to their kinematic 
analogs [g?]. This comparison is important as similar behavior suggests that 
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the dynamical control behavior has its origins at the simple kinematic level. 

The attractive topology of the quantum control landscape, which will be 
affirmed in this work, may be expected to contribute to the generally observed 



favorable lack of scaling of search effort with N [31|, l33| . The attractive global 
topology, however, does not preclude the possibility that complex local land- 
scape structural features may influence the required search effort, particularly 
when using a local search procedure such as a gradient algorithm. The high 
dimensionality of the control landscape (here, the dimensionality is nominally 
infinite as s{t) is a continuous function) renders the direct study of its local 
structure difficult, but useful information about the local landscape features 
can be obtained by examining the trajectories taken during a search from an 
initial to final control. Ultimately, the goal is to understand how the underlying 
control landscape determines the scaling of the required search effort with N. 

The remainder of this work is organized as follows. Section [5] formulates the 
quantum control problem, defines relevant landscape structure metrics, outlines 
the optimization procedure, and defines the model quantum systems. As a 
baseline reference to the optimizations. Section |3] presents the statistical distri- 
butions of Pi^f values obtained when random control fields are applied. Section 
0] shows the important result that no traps were encountered upon optimization 
of Pi^f in ~5000 test cases. Section [S] presents optimization results over varying 
control targets, Hamiltonians and control fields, with the additional general re- 
sult that the search effort is invariant to the system complexity characterized by 
TV, although the absolute search effort varies widely for different circumstances. 
In Section |6l the effect of landscape features on search effort is explored for the 
optimal searches performed in Section [5] using the metrics defined in Section [2j 
Finally, Section [7] presents concluding remarks. 



2 Methods 

2.1 Formulation of the Control Objective 

Consider a quantum system of N levels . . . , |A^) whose dynamics are driven 
by the time-dependent Hamiltonian H{t) — Hq — fJ,e{t), where Hq describes the 
free dynamics of the system, /i is the dipole operator, and e{t) is the control field. 
The time-evolution of the quantum system is given by \ip{t)) — U{t,0)\ip{0)), 
where U{t, 0) is the unitary evolution matrix covering the dynamics from time 
t — to time t and IV'(O)} is the state of the quantum system at i = 0. The 
dynamics of U are governed by the time-dependent Schrodinger equation 

th^^^^ = H{t)U{t,0), (7(0,0)=!. (1) 

The control objective is to maximize the transition probability P^-j. / of pop- 
ulation transfer from an initial state \i) to a target state |/) of the system at 

tllTLG T' 

P,^fiT)^\{f\UiT,Om\ (2) 
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The variation of Pi^f{T) with functional changes in the Hamihonian H{t) is 
obtained by considering smah responses in the propagator U(t, 0): 

d 

ih—SU{t,0)^H{t)SU{t,0) + SH{t)U{t,0), 5J7(0,0)=0 (3) 
SP^^fiT) = {t\6UHT,0)\f){f\U{T,0)\^) + (z|C/t(r,0)|/)(/|JC/(r,0)K). (4) 

Equation ([3]) can be integrated ^] to give 

dU{t,0) = -l f Uit,t')SH{t')U{t',0)dt', (5) 

and substitution of Eq. ([5]) into Eq. dH) yields 

6P,^f{T) = ^Im^ {i\6UHT,0)\f){f\U{T,0)uHt,0)6Hit)Uit,0)\i)dt. (6) 

Within the dipole formulation, SH{t) — —fiSe{t), which gives the functional 
derivative SPi^f /Se{t) from Eq. (6) as 

= ^Im[(*|[/t(i,OV[/(t,0)C/t(T,0)I/)(/|[/(r,0)|z)]. (7) 

We assume that the system is controllable, such that any arbitrary unitary 
matrix U{T, 0) can be generated by a suitably chosen field e{t) at a sufficiently 
large final time T. This condition is equivalent to the requirement that the Lie 
algebra generated from Hq and fi forms a complete set of operators ^2] and T 
is large enough to avoid hindering the dynamics. In general, we may assume 
controllability of an arbitrary quantum system, as uncontrollable quantum sys- 
tems have been shown to constitute a null set in the space of Hamiltonians 
(s^ . Upon satisfaction of the controllability requirement, analysis of the global 
control landscape topology of Eq. ^ with kinematic variables [sij reveals that 
the landscape has no false extrema; the only critical points occur at perfect 
control, Pi^f — 1, and no control, Pi^f = 0. Upon satisfaction of the Jaco- 
bian SU{T,0) /Se{t) being full-rank, the dynamical landscape also has no traps 
[1, 111] and the desired landscape value Pi^ / = 1 corresponds to a submanifold 
of o ptim al fields, which makes the control solutions robust to fiuctuations in e{t) 
(33l. |34| . The latter property is particularly important for laboratory quantum 
control, as it allows for maintaining good yields despite laboratory noise. In 
practice, the rank of 6U{T,0)/Se{t) may be reduced to some degree with no 
impact on the controlled dynamics, as there can still be many readily traversed 
pathways from \i) to |/). However, traps may arise for so-called singular control 
fields where the above Jacobian is significantly rank-deficient. Such situations 



have been known to occur when e{t) — constant is employed [36|-|38j, but this 
situation is generally not physically relevant in the laboratory. Thus, one goal 
of the simulations in this work is to establish whether traps may be encountered 
in optimizations starting from physically reasonable control fields. 
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2.2 Measuring Landscape Structure 

The global landscape topology summarized above provides important information 
about the feasibihty of achieving optimal control. The claimed lack of traps 
means that a control producing a perfect yield can be found starting from any 
initial search point on the landscape (i.e., a point on the landscape corresponds 
to a particular field and its associated transition probability) using a suitable 
hill-climbing algorithm. The validity of this topology in OCT simulations will 
be assessed in this work. 

The presence of a favorable landscape topology does not preclude the pres- 
ence of increasingly complex landscape features as N rises, which could cause 
an increase in the search effort to find a control that gives perfect yield. Thus, 
an understanding of local landscape features (i.e., non-critical point structures) 
is necessary in order to explain and predict the scaling of search effort with 
system complexity. In this work, the local features of the control landscape 
are codified by specific metrics recorded along the search trajectory followed 
from the initial to optimal control field. On a given search trajectory, we may 
parametrize the field e{t) by an index s > to track the progress to the top of 
the landscape. The field starts out at s=0 with e{0,t) and progresses in steps 
s s + ds (i.e., e{s,t) e(s -I- ds,t)) until the trajectory ends at an optimal 
control, Eopt = e(sM,i) at s = sj\/. 

For the purpose of describing the local landscape features, we define (i) 
a distance metric between two fields e{s,t) and e{s',t) {t e [0,r]) based on 
\\e{s,t)—e{s',t)\\, where ||-|| implies an integration over time, and (ii) structure 
metrics based on a Taylor expansion of Pi^f around a field e(s, t) at points on 
the landscape. Analogous metrics of local landscape features were defined in 



61[ using kinematic control variables (i.e., without reference to the dynamics 
of any particular Hamiltonian) and were found to correlate with the observed 
scaling of the search effort with N. From this experience, these metrics are used 
here to provide information about how the features of the landscape determine 
the required search effort using dynamic variables. 

The complexity, or gnarled character, of a search trajectory in control space 
must take into account both the Euclidian distance between the initial and 
final control fields and the actual path length followed from the initial to final 
control over the course of a search. A metric defining this complexity may 
be characterized by the ratio of the trajectory path length ||Ape(i)|| to the 
Euclidian distance between the initial and final control fields ||ABe(t)||, 

2\ 1/2 



^ _ l|Ap£(t)| 



de{s,t) 



ds 



(^j'^ dt[e{sM.t)-e{Q,t)] 



>1 (8) 



The closer is to unity, then the more direct the path, i.e., the closer the 
path is to a straight line in the space of controls being searched over. Following 
a direct path from the initial to optimal control field should result in efficient 
searching, especially by simple local algorithms, because the search trajectory 
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could avoid taking detours along the way to finding an optimal control field. 
This prediction will be assessed in the simulations. 

The local structure metrics of the landscape provide information about what 
the search algorithm "sees" at a particular point on the landscape and may be 
expressed through a Taylor expansion of the cost functional Pi^f , 



VP,^f{s,t)Seis,t)dt+l [ 
^ Jo 



H{t, t')5e{s, t)6e{s, t')dtdt'+- 
(9) 



where VPi^f{s,t) = SPi^f /Se{s,t) is the gradient vector. The structure met- 
rics will be extracted from the kernels of the integrals in Eq. [9l Each metric will 
be labelled by m to indicate its evaluation at the point Sm on the landscape. 
The first-order term in Eq. ^ specifies the slope metric Sm, 



VP,^fiSm,t) 




(10) 



The slope metric is equivalent to the magnitude of the gradient on the landscape 
at the point Sm- Intuitively, a greater value of Sm should result in a locally faster 
ascent due to a more rapid improvement of the yield when taking a step in the 
direction of the gradient. Thus, it is expected that the slope metric may be 
correlated to the observed search effort. 

Additional information about local landscape features can be gained by ex- 
amining the second-order term of the Taylor expansion in Eq. ([9]) , or the Hessian 
matrix, whose elements labelled by t and t' are [s^ 



n{t,t') 



5e{t)5e{t') 



=2Re[(z|f/(0, T)\f) {f\U{T, t)fiU{t, t')fiU{t', 0)\t) 
-{^\UiO,t)^lU{t,T)\f){f\U{T,t')^lU{t',0)\i)l 



t > t'. 
(11) 



The Hessian matrix is symmetric, i.e., Hlt^t') = H^t' ,t). Two simple metrics 
based on the Hessian matrix can provide insight into the landscape structure, 
particularly at the bottom and top of the landscape. The first metric is the 
Hessian trace, 

Tr-H= ( n{t,t)dt, (12) 
Jo 

and the second metric is the curvature of the landscape at a point to, 

2 

/ 1 \ 

r 



VPj^/(s„,t) 



) T T 
dtj^ dt'VP,^fiSm,t)^n{t,t')VP,^f{Sm,t'), 

(13) 

which may be calculated anywhere including near, but not at, the bottom or top 
of the landscape where 'VPi^f{sm,t) = 0. The curvature defined by Eq. is 
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the Hessian projected along the normahzed local gradient direction. Intuitively, 
a larger (positive) value of the Hessian trace and curvature near the bottom of 
the landscape should induce fast climbing [63| . Similarly, a large (negative) value 
of the curvature Cm and Hessian trace TtH near the top should also accelerate 
the approach to the optimum. 



2.3 Optimization procedure 

Many different search algorithms may be used to find an optimal field e{si\i,t) 
maximizing Pi^f. One important goal of this work is to assess whether traps are 
encountered upon climbing the landscape; the existence of traps could preclude 
identification of an optimal control field producing Pi^f 1.0. This landscape 
assessment objective specifically calls for a local (i.e., myopic) search method, 
which will stop climbing at a sub-optimal value of Pi^f if a trap is encountered. 
Global search algorithms (e.g., genetic algorithms) may step over traps, making 
them inappropriate for assessing topology. Additionally, the particular choice 
of search algorithm may significantly influence the absolute effort required to 
find an optimal field; this was found to be the case for optimizing Pi^f using 
kinematic controls ;61j, where gradient, genetic, simplex, and coordinate search 
algorithms were compared. Despite the wide variation in absolute search ef- 
fort with the choice of algorithm, the scaling of the search effort with respect to 
system complexity exhibited the same qualitative trends for all algorithms exam- 
ined. Similarly, in OCT studies from the literature, gradient-based algorithms 
typically converge in - 100 iterations 0, % [13, S, Ei, [Hll, [Hl^, while 



non-gradient simplex and evolutionary searches typically require several hun- 



dred iterations [3|, Il3l . |40| . Importantly, these numbers do not appear strongly 
dependent on N . Considering all of the factors above, a gradient algorithm is 
employed exclusively in this work in order to (a) test the likelihood of encoun- 
tering traps, and (b) seek consistency in exploring optimization effort. 

As the control field e(s, t) depends on the variable s labeling the progression 
of the optimization, the landscape value Pi^f{s) = Pi^f[e{s^t)\ depends on s 
through its functional dependence on e(s,i), < i < T. Thus, the change in 
the landscape value Pi-i./ corresponding to a differential change ds is given by 

dPi^f = (^%^) ds, where 

dP.^f _ f^^^SP^deis^ ^^^^ 



ds Jq Se{s,t) ds 

As the objective is to maximize Pi-^f, we have the demand that '^"^^^ > 0, 
e(s,i) satisfies the differential equation 

de{s,t) _ 6P,- 



ds 5e{s,t)' 



(15) 



where the gradient on the right-hand side is given by Eq. ([T]). Carefully solving 
Eq. (fTS)) coupled to the Schrodinger equation ^ is essential for obtaining re- 
liable landscape climbing results, especially for assessing the presence of traps. 
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The present search algorithm, incorporated into MATLAB [6J| , solves Eq. (fT5|l 
using a fourth order Runge-Kutta integrator with a variable step size to deter- 
mine the control field at the next iteration. Of additional special interest here 
is the required search effort, or the number of algorithmic iterations M required 
to reach the desired Pi^f value, when starting from an initial random control 
field. 



2.4 Design of quantum systems for simulations 

The goals of the simulations are to (a) assess whether traps are encountered in 
carefully performed optimizations and (b) explore general trends in the scaling 
of search effort to find optimal controls in relation to system complexity. For 
a proper assessment of goal (a), as well as in the simulations for (b), no flu- 
ence or other direct constraints are placed on the controls, aside from a fine 
time discretization of the field. Since an infinite variety of structures for Hq 
and fi can arise, a thorough sampling of all physically relevant structures is 
infeasible. Nevertheless, a modest number of variations in Hq, fi, and choice 
of \i) and |/) can capture broad classes of physical phenomena. Increasing N 
while holding the \i) — >■ |/) target transition fixed corresponds to exciting the 
same transition in homologous molecules of increasing size. The circumstance 
of fixing N and the target transition while varying the dipole matrix structure 
corresponds to controlling homologous molecules of similar size with different 
transition couplings. Choosing the target transition as |1) — >■ |7V) and increasing 
TV corresponds to exciting larger molecules to an ever receding highest quantum 
level. In practice, the target \i) — > |/) transition, dipole matrix structure, and 
TV will likely vary simultaneously in the laboratory. The results here should 
both provide diverse test scenarios for the presence of landscape traps as well as 
capture the qualitative search effort scaling trends. Comparisons to the corre- 
sponding laboratory situations will be made at relevant points throughout the 
work. 

For all of the simulations in this work, we consider an A'^-level quantum 
system whose Hamiltonian is expressed in arbitrary dimensionless units. Two 
general choices of nondegenerate, diagonal Hq are employed, corresponding qual- 
itatively to a rigid rotor or an anharmonic oscillator. The energy levels of the 
rigid rotor are given by 

AT-l 

where 7 is a constant. In the results presented here, 7 = 0.25, but varying 7 
was found to have no significant effect on the scaling of search effort with N or 
on the local landscape structure. The energy levels of the anharmonic oscillator 
are 

N-l 



'v 



3 + 



(17) 



9 



where uj — 5 and V = 1200 for all results presented here. Variation of uj and V 
were found not to affect the search effort scaling, provided that they were chosen 
to allow for significantly more bound states than the value of N employed in 
the simulations. The above choices of uj and V provide 120 bound states. The 
Ho structures given in Eqs. ((TC)) and PT)) will be referred to respectively as the 
rotor and oscillator Hq structures later. 

Two physically relevant dipole real matrix structures will be considered. For 
many physical systems the coupling between states generally decreases as the 
difference between the quantum numbers of the states increases, and the present 
choices of /i take this property into account. We first choose fi to have the simple 
structure 



(18) 



where Z? G [0,1] is the drop-off rate and all elements of fj, have a random phase 
of ±1 with the restriction that remains symmetric. We further specify that 
Hif = 0, thereby eliminating a direct transition from the initial state \i) to 
desired target state |/). 

In order to generalize the structure of from that shown in Eq. (|18p . we 
alternatively chose /i to have the form. 
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(19) 



The successive superdiagonal elements ai, i — 1 . . . N — 1, are each chosen from 
particular uniform random distributions such that ai S [0.8, 1], 02 S [0.7,0.9], 
as e [0.6,0.8], ... ai>io £ [0,0.1]. While preserving symmetry, all nonzero 
elements have a random phase of ±1, and =0. The choice of the dipole 
matrices in Eqs. (fT8|) and (fT9|) . respectively, will be referred to the D and 
a structures later. The freedom inherent in randomly drawing the coupling 
matrices provides a broad family of systems to assess the landscape topology, 
structural features, and search effort scaling behavior. 

In many OCT studies, the initial control field is chosen based on knowledge of 
the physical system. For example, the component spectral frequencies are often 
picked to be resonant with certain transitions in Hq, or a spectral bandwidth is 
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chosen that encompasses the desired transitions. In this work, the initial electric 
field £(0, t) is discretized on a time interval t € [0, 28] into 2048 time-points. The 
choice of r=28 and 2048 discretized time-points was found to be sufficient to 
resolve the fastest modulation in the field e{s,t) and the fastest modulation in 
the wavefunction \ilj{t)) for all systems of iV < 30. For simulations involving the 
|1) — >■ I A'') transition for TV > 30, 4096 time points were used to ensure sufficient 
resolution. 

The initial field at s = is chosen as 

K 

J2 sin (uJkt + te[0,T] (20) 

k=l 

where (3 is an envelope parameter (in all simulations, /3=0.05), K is the number 
of frequency components, (j)k is a random phase on [0, 2it], and F is the square 
root of the field fiuence. Prior to multiplication by F, the field is normalized to 
have unit fiuence. The frequencies {wfc} are chosen randomly on a pre-defined 
bandwidth with maximal frequency fi. In most simulations, fl corresponds 
to the frequency of the |1) — >■ |/) transition in Hq, but in Section [Ol other 
choices of are employed. Following selection of the initial frequencies {oJk} 
and the field fiuence F, the electric field is allowed to vary freely over the 
optimization in terms of each of its time-points e{s,tj), j ~ 1,2,..., 2048 (or 
e{s,tj), j — 1,2, .. . ,4096 for some cases where N > 30) as control variables 
starting at s = and iteratively moving ahead as s — > s -|- As. 



e{0,t) = i^exp 



3 Statistical Distribution of Pi^j Yields 

It is instructive to examine the statistical distribution of Pi^f values upon 
making random choices for the initial control field e{0,t) because many OCE 
searches for effective controls start with a random trial choice. Of particular 
interest is whether the optimization searches, on average, start at more or less 
favorable landscape values as N increases. 

A detailed mathematical analysis of the Pi^f objective with kinematic con- 



trols shows that the statistics satisfy a /3-distribution [65|- As N increases, this 
distribution becomes skewed towards smaller P^-^f values. This qualitative be- 
havior has also been observed for initial choices of random control fields e{0,t) 



34l | for the target transition |1) \N). Therefore, simply considering the sta- 
tistical distribution for random trials suggests that increasing search difficulty 
may be encountered as N grows. In order to systematically test the validity of 
this conjecture under different initial conditions, we chose (a) target transitions 
|1) -■ |5), |1) -> |10), and |1) \N), (b) control field fiuence i^=10, 1, 0.1 and 
(c) dipole matrices of structure D in Eq. (ITSl) with D=0.5, 0.2, for N ranging 
from 5 to 40. The statistics were obtained for 10^ different randomly generated 
control fields for each set of parameters |1) \ f), D, and F. All control fields 
had if =20 frequencies randomly distributed on the bandwidth with maximal 
frequency fl=LUf, where denotes the frequency corresponding to the |1) — >■ |/) 
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transition. Results using the rotor Hamiltonian given in Eq. (jl6l) are shown 
here; choice of the oscillator Hamiltonian in Eq. ([TT]) produced qualitatively 
similar results. 

Figure [T] presents the distribution functions for the |1) — |A^) transition with 
fields of F=10 for A^=10, 15, and 20, revealing a shift towards reduced values 
of Pi^f as N rises. The inset of Figure [1] shows the mean of the statistical 
distribution versus N for the cases of different targets, field strengths, and dipole 
matrix drop off rates as labeled in the legend, where D and F are denoted for 
each Pi~if target. For any fixed target transition (e.g., |1) — |5)), the mean 
of each distribution is independent of N, and the distributions for these cases 
are indistinguishable as N is varied (not shown). For the |1) |A^) transition, 
the mean Pi-^f value decreases rapidly with rising N (note the log scale), in 
accordance with [66| . indicating that it becomes increasingly difficult to find a 
decent initial yield as N rises for the receding target \N). The average initial 
yield for systems with a fixed target transition, however, should not change 
dramatically as system dimension rises. Instead, the fluence of the initial control 
field and dipole coupling strength appear to determine the initial yield, with the 
trends following intuitive insights. As expected, stronger fields result in a greater 
yield than weaker fields; however, at very strong fields (not shown), this trend 
can reverse due to amplitude spreading over all the states. Similarly, lower 
yields are obtained for systems with weaker coupling indicated by smaller D 
values. 



4 Testing for the Presence of Traps on the Land- 
scape 

Of primary importance for the utility of quantum OCT and OCE is the question 
of whether all searches starting from a random initial field e(0, t) can even find an 
optimal field achieving Pi-^f ^^1 without getting trapped at a suboptimal Pi^f 
value. Under reasonable assumptions, the topology of the control landscape has 
been theoretically shown to contain no suboptimal extrema when the system 
is controllable, no constraints are placed on the controls, and the Jacobian 
SU{T,0)/Se{t) is full-rank d, [si, 3J, 35|. Affirming this attractive topological 



prediction is very important, as special instances of traps can be found j36l-l3£ 
under unusual conditions. For the Pi-±f objective, the OCT literature regularly 
reports excellent results [iHl, with maximum yields of Pi-^f 



0.9 or greater. These results are not definitive for fully testing the landscape 
theory, as fluence or other field constraints are typically present, and special 
computational care may be required to eliminate artificial traps due to numerical 
aberrations. The present calculations paid due attention to all such details to 
provide a large-scale test of the landscape topology predictions for Pi~^f. As 
pointed out in Section [^31 a gradient-based algorithm was used because a local 
search will stop if a trap is encountered. It is important to execute the gradient 
algorithm in a stable fashion for this purpose, so a fourth-order Runge-Kutta 
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procedure was employed. 

This work provides broad systematic evidence that optimization searches can 
achieve a high yield of Pi^f > 0. 999 without encountering suboptimal extrema. 
A total of ~5000 individual optimal searches were performed with a wide va- 
riety of control parameters chosen (c.f., Section for N ranging from 5 to 
100. In order to ensure that no false traps resulted from choices of simulation 
parameters, the control field was allowed to have as much fluence as necessary 
and the final time T was chosen to be sufficiently large so as not to impose a 
constraint. The importance of paying proper attention to all numerical details 
was evident for some of the difficult cases with the target transition of |1) — lA'') 
(see Section [5751 for further details) when N—30 and 40 with the £'=0.5 dipole, 
the rotor Hq, and employing 2048 time-points to discretize the control field. 
Out of the ~5000 tests, 12 of the latter category were "trapped" at yields of 
0.997 — 0.998. However, upon interpolation of the trapped control fields on 4096 
time-points and continued ascent with the gradient algorithm, the demanded 
criterion of Pi^f >0.999 was achieved in these cases. Similar results were ob- 
served for optimization of the control objective of generating a target unitary 
transformation U{T, 0) with a control field e{t) to match some target unitary 
matrix W . This objective may be measured by considering the fidelity function 
J = \ \W — U{T, 0) 1 p. In the latter study, 20,000 tests were performed on quan- 
tum systems with 2-16 energy levels; upon choice of a sufficiently fine time-mesh 



and large T, each optimization converg ed to a fidelity value of J < lO"*^ [67|. 

Collectively, these results indicate that the likelihood of finding traps on 
quantum control landscapes is vanishingly small when starting with reasonable 
control fields, allowing access to sufficiently fiexible controls, and pa ying atten- 
tion to numerical details. This result suggests that the traps in 36l438l| are at 
most an extremely rare occurrence on the landscape, and possibly a null set. 
Another consideration is that many practical OCT and OCE studies may be 
considered as quite successful upon even reaching moderate yields when oper- 
ating with various constraints. Importantly, the landscape principles affirmed 
by the tests here imply that under such conditions the enhancement of control 
resources can open up even higher yields. 



5 Search Effort and System Complexity 

The scaling of the required search effort with system complexity can determine 
the feasibility of performing quantum control on polyatomic molecules or sim- 
ilarly complex systems. Intuitively, the expectation is that finding a suitable 
control field would become more difficult as the size of the system increases, be- 
cause additional control pathways involving a larger number of quantum states 
become accessible. The collective OCT literature, however, suggests that the 
required search effort to find an optimal control is generally on the order of 
~ 10^ iterations, 0, H, H, [13, 15, 17, 3§-5^, and systematic optimization of 



Pi^f using kinematic control variables indicates that the search effort scales 
at most very slowly with A^ Successful OCE studies ranging from control 
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of atoms 13, 25 1 to complex protein molecules 28, 3^ further suggest a prac- 
tical level of invariance of search effort to system complexity. Based on these 
collective findings, we performed optimization of Pi^f on a broad sampling of 
systems ranging from A^=5 to iV=100 in order to determine whether scaling 
invariance to N can be demonstrated systematically using dynamical control 
variables. The effects of changing the dipole coupling strength, the control field 
parameters, and the \i) \f) target transition on the search effort and its 
scaling with N are examined here. 

5.1 Varying Dipole Coupling Strength 

Optimizations were performed for systems with N ranging from 5 to 40 as 
well as 7V=100 for the target transitions |1) — ?► |5) and |1) — |10). Dipole 
structures of D=1.0, 0.5, and 0.2 as well as the a structure were examined, 
with Hq given by Eq. (IT6l) or Eq. (IT7|) . For all simulations, the initial control 
fields of F=l had K=20 frequencies randomly chosen on a bandwidth with 
maximal frequency corresponding to the |1} — |/) transition in Hq. Optimal 
searches beginning from 20 such initial fields were performed for each choice 
of N and dipole structure, with the exception of iV=100, where 10 optimal 
searches were performed. In order to normalize the reported search effort with 
respect to the initial Pi^f yields obtained, the counting of iterations was begun 
at Pi^f >0.001, regardless of the initial yield, and random fields producing 
Pi^f >0.01 were discarded. 

Figure m shows the mean search effort versus N for rotor Hq (Eq. ([TE)) . 
(a)) and oscillator Hq (Eq. (H?]), (b)) with D^l.O, 0.5, 0.2, and a dipoles 
and the transitions |1) — |5) (solid symbols) and |1) — > |10) (open symbols). 
Representative statistical error bars are presented for one value of N for each 
choice of D and / . Error bars for other N (with the exception of the smallest 
TV for the oscillator Hq structure) were of similar magnitude. Examination of 
Figure [5] shows two striking trends. First, the search effort for any choice of 
dipole structure is invariant to N, at least for TV > 10 for the |1) — >■ |5) transition 
and > 15 for the |1) |10) transition. This result agrees with earlier work 



using kinematic control variables (6l|. Second, for the same dipole structure 
and target transition, the oscillator Hq structure requires a greater effort than 
for the corresponding conditions with the rotor Hq when the dipole coupling is 
weak {D < 0.5). This result shows how the choice of Hq produces landscapes 
with different local structures, as will be reported in Section |6l 

A more detailed examination of Figure [2] reveals two further trends. Stronger 
coupling (i.e., D^l and a dipoles) results in more efficient searches. This in- 
tuitive result can be explained in terms of the accessible mechanistic pathways 
connecting \i) and |/). With strong coupling, both "ladder climbing" (i.e., 
transitions between adjacent states) and quasi-direct transitions are accessible, 
making it easier to find an optimal field that exploits one of many pathways 
from \i) to \f). With weak coupling, accessibility of only adjacent transitions 
limits the number of pathways, thus making it more difficult to find a field that 
utilizes one of them. This phenomenon is illustrated in Figure [31 which shows 
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the population of each state |1) through |10) of a 10-level system plotted ver- 
sus time, with the goal to transfer all population to |10) at T—28. In Figure 
Elja) (Z)=0.2), each intermediate state |2) through |9} is accessed sequentially 
in going from |1} — J> |10). All such plots for D=0.2 showed involvement of each 
intermediate state. In Figure[Hb) (-D=1.0), only states |2) and |8} are involved; 
the remaining intermediate states were never populated more than 10%. Other 
plots for D=1.0 showed between one and eight intermediate states involved, 
indicating more accessible pathways between [1) and |10). Finally, for both Hq 
structures, the more distant |1) |10) transition generally requires more effort 
than the closer |1) |5) transition, except when D=1.0, where the effort is 
similar. This result can be understood in terms of the dipole coupling as well. 
When 13=1.0, all transitions are equally allowed, so changing the final target 
state does not affect the number of accessible mechanistic pathways, resulting in 
no increase of search effort. Weaker coupling, however, closes off pathways be- 
tween non- adjacent states, further reducing the number of accessible pathways 
as the distance between the initial and final states is increased. Plots of state 
population versus time similar to Figure |3] confirm this behavior (not shown). 

5.2 Varying the Initial Control Field 

In order to isolate the effects of varying the initial control field on the required 
search effort, the rotor Hq (Eq. (IT51) ). a dipole structure and the |1) |5) tran- 
sition were fixed. Sets of 20 simulations were performed for initial field strength 
-F=0.1, 10, and 100 with a bandwidth bounded by D,=ll)^ in order to determine 
the effect of initial fluence (i.e., F^) on search effort; the fluence was allowed to 
vary freely during the landscape ascent. Fields of initial strength F=l with a 
fixed maximal frequency Q—UJ20 as well as an TV-dependent maximal frequency 
Q=ljn/2 were also chosen in order to determine the effects of providing more 
bandwidth than necessary. 

Figure [3] presents the mean search effort versus N with representative sta- 
tistical variation bars. The effort is similar for i^=0.1 and F—1 (included as 
a reference), and i^=10; these searches are the most efficient. Further increas- 
ing the field strength leads to greater search effort: at F=100, the effort scales 
exponentially for TV > 10 (note the least squares line and the log scale of the 
ordinate). This result appears to arise because a strong field can easily spread 
an initial amplitude out among many states, making it difficult to then gather 
all of the amplitude into the target state |/). This conclusion can be verified by 
examining the matrix with elements {|C/j/(T, 0)^} produced by the initial and 
optimal electric fields. When F—1, the initial matrix {\Uif{T, 0)^} is nearly di- 
agonal, since the far off-diagonal elements (including the desired (5,1) element) 
are close to zero, as shown in Figure [5Ka). In contrast, when F=100 in Figure 
[S^b), the initial matrix {|C/j/(T, 0)^} contains many significant off-diagonal el- 
ements, indicating that the amplitude is spread out through many states. The 
{\Uif{T, 0)p} matrices produced by the optimal fields retain the predominantly 
diagonal structure for F—1 and the significant off-diagonal elements for F=100 
(Figure [51 bottom row). When the bandwidth provided is more than necessary 
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to make the |1) |5) transition, the effort grows very slowly with N. The slight 
increase in effort compared to using the maximal frequency fl=Lj5 suggests that 
additional access to unneeded ancillary states makes it more difficult to gather 
all of the amplitude in the target final state; examination of the {[[/^/(T, 0)|^} 
matrices for these cases verified this behavior (not shown). 

5.3 The |1) ^ |A^) Transition 

The simulations above employed a fixed choice of \i) and |/) as N was increased. 
Specifying |1) — > |iV) as the target transition causes the final state to recede from 
the initial state as TV is increased. To accommodate the increasing demands of 
transferring amplitude between successively more distant states, the strength 
of the initial fields was chosen as F=10 and the frequencies were chosen on an 
TV dependent bandwidth with maximal frequency fl=ujN/2. Because the initial 
population in \N) drastically decreases with rising N (c.f.. Figure [Ij, iterations 
were counted starting when the yield reached Pi_>./=0.001 to normalize the 
effort against this discrepancy. 

The results of simulations using D—1.0, 0.5 and a dipole structures with both 
rotor and oscillator Hq structures are shown in Figure [SI The scaling behavior 
with N changes significantly depending on the dipole structure. When D—1.0, 
the effort is invariant to N. Although the distance between the initial and final 
states is rising with iV, when all transitions are equally allowed, the number 
of possible pathways between |1) and \N) is large enough to permit efficient 
optimization even at large N. In contrast, for the a and D=0.5 dipoles, the 
effort scales exponentially with N, as shown with the least squares fit lines on 
the semi- log plot in Figure |6l The 12 falsely trapped cases mentioned in Section 
Hlwere for these simulations employing 2048 time-points with D—0.5 and A^=30 
and 40. The additional resolution gained upon interpolation of the control field 
on 4096 time-points eliminated these false traps with further climbing iterations. 
These iterations were added for computation of the mean search effort in Figure 
[6l Receding target objectives with increasing system complexity (i.e., illustrated 
here with |1) — > \N)) are generally not the case for laboratory OCE, thereby 
evidently avoiding the exponential scaling of effort. 

The observed systematic invariance of search effort with respect to N over 
a wide range of Hamiltonian and initial control field structures verifies that the 
search effort for population transfer does not depend on the system complexity, 
as was the case for kinematic controls [6l|. This result is valid upon making 
a rational choice of the control objective and initial field (i.e., for fixed target 
transition and reasonable initial field strength). The results suggest that under 
such circumstances, controlling complex quantum systems with many degrees of 
freedom should be no more difficult than controlling simple systems. Evidently 
the same conclusion applies to performing OCE for various objectives, where 
the search effort appears to be essentially the same regardless of the system 
complexity when operating with physically appropriate controls 0X i25|, i28|, . 
The next section will address the relationship between the observed trends in 
search effort and the underlying control landscape structure. 
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6 Search Effort and Landscape Structure 



Examination of the relationship between the structure of the control landscape 
and the required search effort makes it possible to obtain further insight into 
the scaling results obtained in Section [5l In this section, we determine the local 
landscape structure in terms of the metrics defined in Section 12.21 Here, the 
notion of structure refers to landscape features other than topological critical 
points; the landscape theory predicts critical points only at Pi^/=0 and 1, 
which was verified by the observed lack of traps in Section |4l 

6.1 Search Trajectories on the Control Landscape 

We first consider the relationship between the search effort and the complexity 
of the trajectories over the landscape taken during the optimal searches using 
the ratio metric i?e defined in Eq. ([8]). The mean values of were calculated 
for all the searches performed in Section [5l Select examples with the rotor Hq 
structure (Eq. (fT6|) ) are plotted in Figure [T] When the search effort is invariant 
to N (i.e., the |1) — > |5) transition with F—1 for a and D=0.2 dipoles in Figure 
[7l), the ratio Re is also invariant to N, in agreement with kinematic results 
|6ll |. In contrast, when effort increases with N (e.g., the |1) — >■ |iV) transition 
or large strength F), the path length correspondingly rises with N. For all 
conditions where search effort is invariant to iV, the ratio R^ is correlated to 
the search effort, as shown in Table 1 for simulations using the rotor Hq (left 
of double line) and oscillator Hq (right of double line); ratios are significantly 
higher for the oscillator Hq, although these do not scale with N. The values of 
the distances ||A£;£|| and ||Ape|| used to define R^ follow the same correlations 
with effort. The differences in values of Re between optimizations using the 
rotor and oscillator Hq structures for weakly coupled dipoles can be explained 
by examination of the landscape slope, as discussed below. 

6.2 Landscape Slopes and Search Effort 

The magnitude of the gradient Sm provides valuable information about how 
fast a search algorithm may improve the yield. Intuitively, a steep slope would 
be conducive to efficient optimization because the yield may improve rapidly 
upon taking an algorithmic step, while a very shallow slope should slow the 
optimization. 

For the optimizations in Section O the slope metric Sq at the initial ran- 
dom control field (or at the first iteration where Pi-^f >0.001) and the point 
of maximal slope Smax were recorded; at Pi^f ~0.001, the slope metric Sq is 
typically small. Both the initial Sq and maximal Smax slope metrics along an 
optimization may be expected to correlate with the required search effort. Fig- 
ure [8] shows the mean value of the initial slope metric Sq (filled symbols) and 
maximal slope metric Smax (open symbols) for selected optimizations from Sec- 
tion[5] The initial and final slope metrics are independent of N under conditions 
where the search effort is also invariant, while both metrics for the |1) — J" \N) 



17 



transition decrease as A'' rises, in accordance with the increase in search effort. 
AU conditions where the effort was dependent on N exhibited the behavior of 
decreasing slope metrics as N rises. For the cases invariant to N, more difficult 
optimizations (e.g., optimization with a weak dipole) have smaller initial and 
maximal slope metrics than easier optimizations, as shown in Table 1. Thus, 
the search effort follows the intuitive conjecture that a steeper slope results in 
more efficient optimization, as was found using kinematic control variables [61 1. 
In general, the linkage of search effort to the gradient depends on the choice of 
search algorithm. Most OCT studies use gradient algorithms, so in such cases 
the search effort may be expected to depend on the initial and/or maximal slope 
metric. However, other "smart" algorithms (e.g., with stochastic logic) can also 
exploit the favorable slopes and direct pathways to the optimum with being 
small. 

An exception to the simple search effort correlation with the initial and 
maximal slope metrics arises for searches using weakly coupled dipoles when 
comparing the two Hq structures with otherwise identical search conditions. The 
effort for the oscillator Hq is drastically higher than for the rotor Hq, but the 
initial and maximal slope metrics are of similar magnitude, as shown in Figure 
[5] and Table 1. This discrepancy can be explained by examining the trajectory 
of the slope metric and the ratio over the course of an optimization. As 
an example, these trajectories for searches with A^=20, D=0.2 and |1) — !> |5) 
transition are compared for the two different Hq structures. Figure [HI shows 
the trajectory of the slope metric (a) and the trajectory of Re (b) for two 
searches with each Hq structure. The trajectories of the slope metric S for 
the rotor Hq share the simple structure of starting near zero at the initial field 
with Pi^f ~ 0.001, rising to a maximum around Pi^f ~ 0.5, and decreasing 
towards the optimum. Similarly, the trajectories of R^ for these searches show a 
simple monotonic rise with Pi^f. In contrast, the trajectories of searches using 
the oscillator Hq structure show a more complex behavior over the landscape. 
Instead of reaching a high at Pi^f ^ 0.5, the maximal slope metric for the 
oscillator searches occurs below Pi^f ~ 0.3, and the slope decreases rapidly 
thereafter. Examination of R^ at Pi^j values (b) corresponding to the rapidly 
decreasing slope metric in (a) shows a fast jump in i?e with Pi^f, indicating a 
relatively "gnarled" landscape region. Finally, the slope metric for the oscillator 
searches drops quickly for Pi^f > 0.8, and the ratio Re rises accordingly. Other 
trajectories for searches using the oscillator Hq with a weakly-coupled dipole 
show similar features, suggesting that an oscillator Hq structure with a weakly 
coupled dipole inherently creates a more gnarled landscape than a rotor Hq with 
the same dipole. 



6.3 Second Order Landscape Structure 

Examination of the second-order landscape structure metrics can provide further 
insight into contributions to the relative search effort required under different 
optimization conditions. Calculations of the Hessian matrix and associated 
structure metrics at the bottom and top of the landscape were performed on 
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the rotor Hq structure for (i) the |1) — >■ |5) transition with the a and D=Q.5 
dipole structures for F=l, (ii) the a dipole structure for F=100, and (iii) D=0.5 
for the |1) \N) transition with F=10. With the oscillator Hq structure, the 
calculations were performed for the D=0.5 dipole and F=l with |1} |5) 
transition. In order to obtain Hessian matrices reliably representing the bottom 
and top of the landscape, all optimizations began at Pi^f < 1 x 10~^ and the 
convergence criterion was Pi^f > 0.99999. 

It has been shown theoretically that the Hessian spectrum at the bottom of 
the landscape has at most two nonzero positive eigenvalues and the spectrum at 
the top contains at most 27V- 2 nonzero negative eigenvalues 33j. This analysis 
is verified by our numerical results. Figure [TU] shows the Hessian spectra at the 
top of the landscape for individual optimizations of |1) — )■ |5) transition with 
rotor Hq structure, F=l, and D=1.Q (Figure [TOTa). for N ranging from 5 to 
30) and 15=0.5 (Figure ITUT b) . for N ranging from 5 to 15). The vertical dotted 
lines denote the eigenvalue index of 2A''-2 for each N reported. In the case 
of D=1.0, there is always a clear distinction between the (27V-2)th eigenvalue 
(~ —10) and the (2A^-3)th eigenvalue (> —0.01). The magnitude of the largest 
and smallest nonzero eigenvalues does not change with N. For _D=0.5, the drop 
in eigenvalue magnitude at the index 2A^-2 is apparent at N—5 and 10 (note log 
scale on the ordinate in Figure [TUTb)). By A^=15, the distinction between the 
final nonzero and first zero eigenvalue is expected to occur between the 28th and 
29th eigenvalues, however the eigenvalues are already of very small magnitude 
by the 23rd eigenvalue. Recording the eigenvalues for larger values of N with 
D=0.5 revealed similar patterns of eigenvalue behavior. This result shows that 
for large N with weak dipole couplings, fewer than 2A^-2 negative eigenvalues can 
be expected at the top of the landscape, and there is no clear boundary between 
the zero and nonzero eigenvalues. Fewer than 2iV-2 nonzero eigenvalues were 
also observed for TV > 15 using the oscillator Hq structure with D=0.5 (not 
shown). With strong dipolar couplings (i.e., D=1.0), there are always exactly 
2A''-2 nonzero eigenvalues; for the a dipole structure, exactly 2A^-2 eigenvalues 
persist through iV=30, and by A^=40 there are fewer than 2A^-2 eigenvalues (not 
shown). At the bottom of the landscape, there is a clear distinction between 
the two positive eigenvalues and the remaining zero eigenvalues, which occurred 
under all search conditions (not shown). These observations about the Hessian 
eigenvalues at the bottom and top of the landscape validate the theoretically 
predicted spectra [sl] . Additionally, the number of non-zero Hessian eigenvalues 
at the top of the landscape influences the robustness of the control outcome to 
field noise; the presence of fewer such eigenvalues enhances the robustness [63 1. 

Examination of the Hessian trace and curvature metrics (c.f. Eqs. (jl2p 
and (US])) at the bottom and top of the landscape yielded intuitive correlations 
between these metrics and the required search effort, as was the case with the 
slope metric. As graphs of these metrics versus iV are similar to Figure |S1 
the data are not plotted again. Near the bottom of the landscape, both the 
Hessian trace and curvature metrics are invariant to N when the search effort is 
also invariant, and smaller values of these metrics are recorded for more difficult 
search conditions (e.g., oscillator Hq, small dipole coupling). Where exponential 
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scaling of search effort with N was found, both metrics decrease exponentially 
with N near the bottom of the landscape. At the top of the landscape {Pi^f > 
0.99999), the Hessian trace is proportional to N, regardless of search parameters, 
due to its dependence on the dipole norm ||/z||^ ^333. The curvature exhibits 
intuitive correlation with the search effort, remaining constant with N for cases 
that lack search effort scaling, and decreasing in magnitude with N where scaling 
is observed. Thus, all of the landscape structure metrics examined in this section 
correlate in an intuitive way with the required optimization search effort. These 
results show that the landscape structure metrics provide a good method to 
predict the relative required search effort under a variety of conditions. 



7 Conclusion 

This work addressed two major issues surrounding optimal control of population 
transfer in quantum systems. The first objective explored the fundamental topic 
of whether suboptimal trapping extrema are encountered while searching for an 
optimal control field. The second objective examined how the required effort to 
find an optimal control field scales with the complexity of the quantum system 
as measured by its size TV. 

The possible existence of traps on the control landscape is of both basic and 
practical importance. Quantum control landscapes can rigorously be shown 
to contain no traps under simple physical assumptions [3l|, [ss H sal. The vast 
OCT literature supports the ability to reach excellent yields [l|-|3|, iM Hi- 
Is^, although these works are not definitive with regard to the landscape due 
to control field constraint s ty pically being present. The recent identification 



of trapping conditions |36l - l38| under unusual circumstances necessitates a more 
explicit investigation of whether traps can be expected when performing normal 
optimizations. 

The simulations in this work found no evidence of trapping behavior on the 
control landscape for Pi^f ■ Of the ~5000 searches performed, a total of 12 were 
initially found to be putative traps warranting further investigation. Enhancing 
the time resolution established that the latter traps were in fact false, with 
all optimization searches then reaching Pi^j > 0.999. The identification of 
false trapping behavior due to numerical constraints illustrates the need for 
special care in performing simulations and the general need for due attention 
to all physical constraints on the field dynamics when a high yield is desired. 
The lack of observed traps on the Pi^f landscape is consistent with results 
reported for the landscape corresponding to the generation of arbitrary unitary 
transformations U{T,0), where ^20,00 op timizations were performed, all of 
which reached an optimal fidelity value |67| . 

The second issue studied here of search effort scaling with N is primarily of 
practical importance, indicating whether the control of large, complex quantum 
systems in the laboratory is feasible. The OCT literature collectively suggests 
that the required search effort to find an optimal control may be independent 
of the complexity (i.e., here captured by N) of the target quantum system 
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0, S, i, H El, [H [ai-IHi] . The results from this work systematically verify this 
behavior and identify the control conditions sufficient for the search effort scaling 
to be independent of A^. Specifically, choosing a fixed target transition \i) —>■ \ f) 
results in the scaling of effort being invariant to N across a wide range of dipole 
matrix structures and reasonable initial control field parameters, although the 
absolute search effort can vary widely. This attractive behavior breaks down, 
however, upon choosing targets that themselves increase in complexity with the 
system (e.g., |1) — >■ |A^)) or starting with a large initial control field strength 
for a fixed target transition, where the wavefunction amplitude spreads widely 
before finally being drawn into the target state. 

The observed search effort was found to correlate with the landscape features, 
as measured by the distance and structure metrics. The scaling of the ratio of 
path length to Euclidian distance Re with N follows that of the search effort; Re 
only increases with N for the difficult cases such as the |1) — J> \N) target or with 
a large initial field fluence. For cases with scaling invariant to N, the relative 
search effort can be predicted by the value of Re, with greater values of this 
metric correlating with a greater search effort. Analysis of the local structure 
of the landscape shows that the search effort correlates with the slope metric 
(gradient norm) in an intuitive manner. A steeper landscape slope both at the 
initial control field and at the point of maximal slope results in a lower search 
effort than a shallow slope. The landscape slopes at these points are invariant 
to N, except for the cases where the search effort scales with A^, for which 
both initial and maximal slopes decrease as N rises. A similar correlation of 
search effort with the curvature metric near the bottom and top of the landscape 
with N was observed. Finally, the collective dynamic findings on search effort 
show a strong relation to analogous behavior found using kinematic variables 



6l| . Although clearly additional dynamical features occur (e.g., through the 



amplitude and structure of the dipole couplings), much of the basic invariant 
scaling findings with N appear to have their origins in the underlying simpler 
kinematic control formulation. 

This work addressed many classes of control Hamiltonians in order to demon- 
strate the broad applicability of the two main results in this work. However, 
some classes of quantum systems, such as those containing degenerate energy 
levels or additional symmetry in the dipole matrix, were not addressed here. 
Provided that such systems are controllable [sl] (e.g., where dipole couplings 
break the symmetry produced by degenerate states), the favorable topological 
and scaling results are expected to hold. For other special classes of systems 
that are uncontrollable or nearly so (e.g., a harmonic oscillator), special care in 
the choice of controls may be needed to avoid traps on the landscape arising 
from the lack of system controllability. Most classes of quantum systems, how- 
ever, are expected to satisfy the controllability requirement and thus exhibit 
qualitatively similar behavior in terms of landscape topology and search effort 
seen here. 

The favorable scaling of Pi^f with N suggests that optimization of state 
preparation with a suitable set of controls should be relatively easy to attain 
using OCE, even with complex systems. Although the quantum systems em- 



21 



ployed here do not model any particular real system, the results using the rigid 
rotor and anharmonic oscillator Hq structures indicate that some quantum sys- 
tems may generate a landscape with a more gnarled local structure than others, 
leading to wide variations in the absolute search effort required to find an op- 
timal control. Nevertheless, a family of quantum systems that are difhcult to 
optimize may still show invariant scaling with N. These results are consistent 



with successful OCE studies on complex molecules such as proteins [28|,l30|, even 
though the laboratory conditions are more involved than the ideal circumstances 
presented here. 

Overall, this work demonstrated that both the topology and the local struc- 
ture of the control landscape for population transfer are conducive to efficient 
optimal control. Extensive simulations did not encounter traps on the landscape 
upon reasonable choices of Hamiltonians, initial control fields, and careful nu- 
merical optimization. The invariance of scaling of the search effort with system 
complexity was shown to be due to favorable local landscape structure that does 
not grow more complex with system size N. Besides state preparation, recent 
studies generalize these landscape topology, features, and optimization scaling 



results to the preparation of unitary transformations [67[ and broader classes of 
observablcs 68]. 
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Captions 

Table [TJ Mean search effort, ratio R^, initial slope metric Sq, and maximal 
slope metric iSmax for all simulations that showed invariant scaling effort to N. 
Recorded values are taken from simulations at 7V=20, but for other N the val- 
ues were similar. The values to the left of the double line are from simulations 
using the rotor Hq (Eq. HI])), and the values to the right of the double line are 
from simulations using the oscillator Hg (Eq. ([17])). A comparison of the land- 
scape metrics with the effort shows that the two are correlated. The "easiest" 
optimizations {D—1.0) have the lowest ratio and the highest initial Sq and 
maximal iSmax slope metrics, while the "hardest" optimizations (D=0.2) have 
the highest ratios Rg and lowest initial and maximal slope metrics. The effort 
using the oscillator Hq is always greater than for the rotor Hq, and the metrics 
show corresponding increases. 

Figure 1. Statistical distributions of Pi^f values for N=10, 15, and 20, 
with D=Q.5 and F=10 for the |1) \N) target. The inset depicts the mean 
value of distributions of initial Pi^f values for different dipoles, targets, and 
field parameters. The target transition, dipole drop off rate D and field fluence 
F are denoted as Pi^f, D, F in the legend. The mean initial value decreases 
for the |1) \N) target, but is constant for fixed target transitions. Statistical 
error bars are shown for the |1) \N) transition, and representative error bars 
for the other cases are shown as well. Some points are shifted on the x-axis for 
graphical clarity. 

Figure 2. Required mean search effort versus N for the target transitions 
|1) — >■ |5) (solid shapes) and |1) — > |10) (open shapes) for Hamiltonians with 
dipole structures of D=1.0 (squares), D=0.5 (circles), D—0.2 (down triangles) 
and a (side triangles), with Hq given by Eq. (fTH)) (a) and by Eq. ITT]) (b) Search 
effort is invariant to N in all cases (excepting some cases where the effort for 
the smallest N recorded is significantly lower than for remaining N), but the 
absolute effort is greater for weak coupled dipoles, the |1) |10) transition, and 
oscillator Hq structure. Some points are shifted on the abscissa for graphical 
clarity. 

Figure 3. Population of states versus time for a 10-level system with target 
|1) — >■ |10). (a) D=0.2, and all intermediate states |2) through |9) are accessed 
sequentially on the way from |1) to |10), consistent with a ladder-climbing mech- 
anism, (b) D=1.0, and only states |2) and |8) play a significant role (all other 
intermediate states are never populated more than 10% and arc not shown). 

Figure 4. Required mean search effort versus N for the target transition 
|1) — ?► 1 5) and a dipole structure with varying initial field strength and band- 
width. The strength (solid shapes) or bandwidth (open shapes) is labeled in 
the legend. For low strength and reasonable bandwidth, effort is invariant to 
N. For high fluence (F=100), effort scales exponentially with N, as shown by 
the least squares fit on the semi-log plot. For large bandwidth range ft, effort 
increases through iV=20 and then levels off. Some points are shifted on the 
abscissa for graphical clarity. 

Figure 5. Plots of the absolute value of the matrix elements of the prop- 
agator {|J7i/(T, 0)p} at initial random control fields (top) and optimal fields 
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(bottom) for 7V=10 and the target |1) — J> |5) transition under the conditions 
F=l (a) and -F=100 (b). The (5,1) element is circled in each plot and has 
a value of 1.0 at the optimal fields and a value of close to zero at the initial 
fields. Both the initial and final {\Uif{T, 0)1"^} matrices are nearly diagonal for 
F=l, while at F=100, many off-diagonal elements are non-zero, indicating that 
the population is spread out among many states. Such a {|t/i/(T, 0)|^} matrix 
structure at i^=100 results in a greater search effort because it becomes more 
difficult to gather all of the amplitude in a single final target state. 

Figure 6. Mean search effort versus N for the |1) \N) transition. When 
all transitions are allowed (D=1.0, squares), the effort is invariant to N. When 
the coupling strength decreases with distance between the states (a, triangles, 
and D=Q.5, circles), the effort scales exponentially with A^, as shown by the least 
squares fit lines on the semi-log plot. Results are qualitatively the same for the 
rotor Hamiltonian (filled symbols) and oscillator Hamiltonian (open symbols). 

Figure \7\ Ratio of the control search path length to the Euclidian 
distance versus N for selected optimizations from Section [Sj The ratio is in- 
variant to N for the |1) — J> |5) transition and small strength F; the effort was 
also invariant to N for these cases. Rg increases with N for optimizations with 
iV-dependent effort (i.e., |1) -4 |A^) and |1) -J- |5) with F=100). Regardless of 
these variations, the values of R^ are generally close to 1, indicating that the 
searches follow direct trajectories in the space of controls. 

Figure [8l Initial (filled symbols) slope metric Sq and maximal (open sym- 
bols) slope metric S,nax versus N for selected cases from Section [5j With the 
exception of the |1) — > |A^) transition, the initial and final slope metrics are 
invariant to N, in agreement with the observed scaling behavior. Some points 
are shifted on the abscissa for graphical clarity. 

Figure [9l Trajectory of the slope metric (a) and ratio R^ (b) for two 
searches using the oscillator Hq (solid lines) and rotor Hq (dashed lines). £'=0.2, 
F=l and the target transition is |1) — >■ |5). The trajectories for searches using 
the rotor Hq are less complex than those for searches using the oscillator Hq . 

Figure IIOI Hessian eigenvalues at the top of the landscape plotted versus 
their index. All optimizations used rotor Hq structure and F—1 for the |1) — )> |5) 
transition, (a) optimizations with D~1.0. (b) optimizations with D—0.5] note 
the logarithmic scale. The vertical dotted lines show the value of 2A^-2 for each 
N, and the 2A^-2 rule is obeyed. In each case a few of the zero eigenvalues are 
shown for graphical clarity. 
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